begin

    large_image = load("philip.jpg")

    image = shrink_image(large_image, 7)

end

size(large_image)

size(image)

kernel = Kernel.gaussian((2, 2))

show_colored_kernel(kernel)

sum(kernel)

kernel_sharp = centered([

        -0.5 -1.0 -0.5

        -1.0  7.0 -1.0

        -0.5 -1.0 -0.5])

sharp_image = convolve(image, kernel_sharp)

sum(kernel_sharp)

heatmap_2d_fourier_spectrum(sharp_image)

image

heatmap_2d_fourier_spectrum(image)

conv_image = convolve(image, kernel)

heatmap_2d_fourier_spectrum(conv_image)

kernel_sobel = Kernel.sobel()[2]

show_colored_kernel(kernel_sobel)

sum(kernel_sobel)

3 * Gray.(abs.(convolve(image, kernel_sobel)))

begin

    big_zebra = load("zebra.jpg")

    zebra = shrink_image(big_zebra, 7)

end

plot_1d_fourier_spectrum(zebra)

heatmap_2d_fourier_spectrum(zebra)

zebra_kernel = Kernel.gaussian((2, 2))

conv_zebra = convolve(zebra, zebra_kernel)

plot_1d_fourier_spectrum(conv_zebra)

heatmap_2d_fourier_spectrum(conv_zebra)

begin

    using Statistics

    using Images

    using FFTW

    using Plots

    using DSP

    using ImageFiltering

    using PlutoUI

end

function show_colored_kernel(kernel)

    to_rgb(x) = RGB(max(-x, 0), max(x, 0), 0)

    to_rgb.(kernel) / maximum(abs.(kernel))

end

function decimate(arr, ratio=5)

    return arr[1:ratio:end, 1:ratio:end]

end

function shrink_image(image, ratio=5)

    (height, width) = size(image)

    new_height = height ÷ ratio - 1

    new_width = width ÷ ratio - 1

    list = [

        mean(image[

            ratio * i:ratio * (i + 1),

            ratio * j:ratio * (j + 1),

        ])

        for j in 1:new_width

        for i in 1:new_height

    ]

    reshape(list, new_height, new_width)

end

function rgb_to_float(color)

    return mean([color.r, color.g, color.b])

end

function fourier_spectrum_magnitudes(img)

    grey_values = rgb_to_float.(img)

    spectrum = fftshift(fft(grey_values))

    return abs.(spectrum)

end

function plot_1d_fourier_spectrum(img, dims=1)

    spectrum = fourier_spectrum_magnitudes(img)

    plot(centered(mean(spectrum, dims=1)[1:end]))

end

function heatmap_2d_fourier_spectrum(img)

    heatmap(log.(fourier_spectrum_magnitudes(img)))

end

function clamp_at_boundary(M, i, j)

    return M[

        clamp(i, 1, size(M, 1)),

        clamp(j, 1, size(M, 2)),    

    ]

end

function rolloff_boundary(M, i, j)

    if (1 ≤ i ≤ size(M, 1)) && (1 ≤ j ≤ size(M, 2))

        return M[i, j]

    else

        return 0 * M[1, 1]

    end

end

function convolve(M, kernel, M_index_func=clamp_at_boundary)

    height = size(kernel, 1)

    width = size(kernel, 2)

    half_height = height ÷ 2

    half_width = width ÷ 2

    new_image = similar(M)

    # (i, j) loop over the original image

    @inbounds for i in 1:size(M, 1)

        for j in 1:size(M, 2)

            # (k, l) loop over the neighbouring pixels

            new_image[i, j] = sum([

                kernel[k, l] * M_index_func(M, i - k, j - l)

                for k in -half_height:-half_height + height - 1

                for l in -half_width:-half_width + width - 1

            ])

        end

    end

    return new_image

end

box_blur(n) = centered(ones(n, n) ./ (n^2))

function gauss_blur(n, sigma=0.25)

    kern = gaussian((n, n), sigma)

    return kern / sum(kern)

end